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In Supplemental Material, we define and explain local 
helicity modulus and local superfluid density. The con¬ 
tents mainly focus on the numerical techniques of quan¬ 
tum Monte Carlo (QMC) simulation. Readers who are 
interested in only physical results may skip this Supple¬ 
ment Material. 

The helicity modulus Ytot and the superfluid density 
^HOMO dimensional uniform system [l| are defined 
by using a winding number vector W 0 as 

where t is the ’’averaged” short-range hopping amplitude, 
L is a linear size of the system. The a-component Wa 
of W stands for the number of winding particles along 
the a axis in the path-integral formalism, and can be 
evaluated by QMC method. The value Itot is usually in¬ 
sensitive to details of the system, but may depend 

strongly on them such as values of hopping amplitudes 
and interactions. In fact, for instance, it is well known 
that the helicity modulus exhibits a universal jump at a 
Kosterlitz-Thouless (KT) transition point in two dimen¬ 
sions. 

Now, keeping in mind the above facts of the uniform 
systems, let us discuss local superfluidity and local helic¬ 
ity modulus in spatially inhomogeneous systems. To see 
the universal nature of the systems, we will mainly use 
the local helicity modulus rather than local superfluid 
density. In order to theoretically define helicity modu¬ 
lus, we have to impose a periodic boundary condition 
along at least one spatial direction, because the modulus 
is defined by the energy variation after slightly twisting 
the boundary condition. Namely, by definition, it is im¬ 
possible to define ’’truly” local helicity modulus at any 
position X on the d dimensional space. However, when 
d > 2, we can introduce a local helicity modulus defined 
at any position on one special x direction by imposing 
a periodic boundary condition along a direction perpen¬ 
dicular to the X axis. For a d-dimensional system with 
d > 2, let us dehne the local helicity modulus Y{x) on 
X axis such that the total helicity modulus is given by 
summation of Y(x) along the x direction: 

= ( 2 ) 


Following Eqs. (HD and m, we then define the relation 
between the local superfulid density ps (x) and Y {x) as 

Ps{x) = (3) 

In this equation, Ps(x) stands for the superfluid density 
along a direction perpendicular to the x axis. For exam¬ 
ple, if we consider a two dimensional system with open 
boundary condition along the x axis, Psix) is the super- 
fluid density along the y direction at x. 

Next we explain how the local helicity modulus is re¬ 
lated to numerical evaluated quantities such as W. Sim¬ 
ilarly to local helicity modulus, winding (topological) 
numbers including W also has a non-local nature and 
a periodic boundary condition is necessary along at least 
one direction to define them. Paying attention to these 
properties, we can introduce the partial winding number 
vector w{x) at each position on the x axis. The summa¬ 
tion of w{x) over the x axis is equivalent to the original 
winding-number vector W: 

W = '^w{x). (4) 


The X component Wx{x) of w(x), which can be evalu¬ 
ated by QMC, does not satisfy the property of winding 
number, since information about the whole range of the 
X axis is necessary to define winding number along the 
X axis. However, sum of them, i.e., Wx can be a wind¬ 
ing number if we impose a periodic boundary condition 
along the x direction. Each remaining component Wa (x) 
{a ^ x) is the number of winding particles along the a 
direction at x, and thus it has a topological nature. 

From Eqs. (HD, (ED and O, we find that w{x) and Ps{x) 
should satisfy the following relation, 


Y{x) = 


{W ■ w{x)) 


( 5 ) 


When an open boundary condition is applied along the 
X axis as in our setup in the main text, the x component 
of the winding vector becomes zero and Ps (x) represents 
the superfluid density along the y axis at the position x. 
In this way, we can evaluate the local helicity modulus 
and local superfluid density (see Eig.2 of the main text). 
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We study the boundary nature of trapped bosonic Mott insulators in optical square lattices, by 
performing quantum Monte Carlo simulation. We show that a finite superfluid density generally 
emerges in the incommensurate-filling (IC) boundary region around the bulk Mott state, irrespec¬ 
tively of the width of the IC region. Both off-diagonal and density correlation functions in the IC 
boundary region exhibit a nearly power-law decay. The power-law behavior and superfluidity are 
well developed below a characteristic temperature. These results indicate that a gapless boundary 
mode always emerges in any atomic Mott insulators on optical lattices. This further implies that 
if we consider a topological insulating state in Bose or Fermi atomic systems, its boundary pos¬ 
sesses at least two gapless modes (or coupled modes) of an above IC edge state and the intrinsic 
topologically-protected edge state. 

PACS numbers: 05.30.Jp,37.10.Jk,73.43.-f 


Introduction.— Topological insulators (TIs) 0,i , more 
widely, symmetry-protected-topological (SPT) states [§- 

, have been vividly studied as new quantum many-body 
states in the last decade. These gapful states cannot be 
characterized by any local order parameter, while they 
generally possess a gapless edge/surface mode. Each 
SPT phase is protected by certain symmetries, namely, 
it is stable against any perturbation keeping the symme¬ 
tries. A complete classification of TIs and the relation¬ 
ship between their bulk symmetry and the corresponding 
surface/edge state have been established for free fermion 
models [fl-Q . Several TI materials have been synthesized 
and their surface/edge states have been observed Si- 

Many physicists stimulated by the study of fermionic 
TIs have been exploring SPT states in spin and bo¬ 
son systems. The Haldane-gap state i [13 of one¬ 
dimensional (ID) spin-1 antiferromagnets is a typical 
SPT state in quantum spin systems, and it is indeed re¬ 
alized in several quasi-ID magnets [ll|, In addition 
to the Haldane state, several SPT phases in ID fermion, 
boson and spin systems have been discussed. In fact, a 
way to classify ID bosonic SPT phases has been proposed 
by tensor product representation [l^ . 

On the other hand, two- or three-dimensional (2D or 
3D) bosonic SPT phases and their edge/surface states 
have been little understood. Several theorists discussed 
the possibility of higher-dimensional bosonic TIs @,0- 
lisjl - and proposed ways to classify them: bosonic TIs in 
spatial dimensions d can be distinguished by a technique 
based on {d -\- l)th group cohomology |15l4l8l| . In those 
studies, some models for bosonic TIs were predicted, but 
it is difficult to realize them in real materials because the 
corresponding Hamiltonians contain various tuned cou¬ 


pling constants. 


For the realization of 2D or 3D SPT phases in boson or 
spin systems, strong interactions among bosonic particles 
or spins are generally necessary. The interaction usually 
makes it quite difficult to analyze the systems, and this 
is a main reason why the theory for 2D or 3D bosonic 
TIs has not been developed in comparison with fermionic 
TIs. Because of the same reason, even non-topological 
(i.e., trivial) gapped phases and their boundary nature 
have not been understood well in the strongly interacting 
boson and spin systems. In a sense, boundary nature of 
gapped states is more important than classification of 
ground states because physical phenomena at boundary 
can be observed and their information often provides a 
experimental way to characterize the bulk state. 


Recently, we have studied a edge state of 2D spin- 


Peierls states IH 
culations Mm 


by quantum Monte Carlo (QMC) cal- 
The Peierls state is a typical trivial 
gapped state in quantum spin systems and it does not ac¬ 
company any spontaneous breaking of basic symmetries 
(such as spin rotation and time-reversal symmetries). We 
showed that if we prepare a sufficiently clean edge of the 
Peierls state with a large enough length (~ 50 sites), we 
can observe a gap less Tomonaga-Luttinger-liquid (TLL) 
like behavior [23| along the edge and the edge spin-spin 
correlation function decays in an almost algebraic fash¬ 
ion. We proposed some experimental ways of detecting 
these gapless edge excitations. In this paper, we will ex¬ 
plore the fundamental nature of boundary states of 2D 
Bose Mott insulator on optical lattices, by QMC compu¬ 
tations. Similarly to the spin-Peierls state and fermionic 
TIs, no spontaneous symmetry breaking occurs and a fi¬ 
nite bulk excitation gap exists in Bose Mott states. They 
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FIG. 1: (color online) (a) Full region of a 2D trapped Bose sys¬ 
tem on a square lattice, (b) Quasi-ID geometry we consider 
and (c) the corresponding confinement potential. Density of 
color stands for the depth of chemical potential. 


have been already realized in ultracold-atomic systems on 
optical lattice 24l-l28l|. Therefore, their boundary nature 
could be an important research subject as a good com¬ 
parison in that of bosonic TIs. 

An important feature of trapped ultracold-atom sys¬ 
tems is that their boundary is always clean. This con¬ 
siderably contrasts with solid systems, whose boundary 
is usually dirty. As a result, we always observe a clean 
and homogeneous boundary region in Bose Mott states. 
We will numerically clarify the boundary properties of 
2D Bose Mott states. They could be experimentally de¬ 
tected in principle. Our findings are also useful to deeply 
understand boundary states of cold-atom TIs as well as 
those of Bose Mott states. 

Model— In this paper, we focus on the 2D soft-core 
Bose Hubbard model with confinement potentials. To 
discuss the finite-size effects systematically, we consider 
systems on the quasi-lD geometry shown in Fig. [TJb). 
This geometry would be hard to be realized in optical 
lattices, but it could be regarded as a boundary part of 
circular or elliptic shaped trapped systems [see Fig.[Tl[a)]. 
The Hamiltonian for the quasi-ID geometry is given as 

x,y 

+U^nly-^V{x)nx,y, ( 1 ) 


where V{x) = /tq + ax^, nx,y = b^ yb^^y and bx,y (bl^y) 
is a boson creation (annihilation) operator at position 
{x,y). Parameters t, U, and V(x) denote hopping am¬ 
plitude, on-site repulsion, and axial confinement poten¬ 
tial along the x-axis direction, respectively. To realize 
the quasi-lD geometry, we impose the periodic (open) 
boundary condition along the y {x) direction. In our 
computations, we fix D/t = 20, at which the bulk system 
can belong to the Mott-insulating state with a single bo¬ 
son per site [ii-ii. We also fix the x-direction length 
Lx = 48, but tune the y-direction length Ly. Increase 
of a means the growth of potential slope between bulk 
Mott and vacuum (empty) regions. 
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FIG. 2: (color online) (a)-(c) Density and (d)-(f) stiffness (oc 
superfluid density) profiles estimated by QMC calculations at 
extremely low temperatures T. Pairs [(a),(d)], [(b),(e)] and 
[(c),(f)] are the results of CASEs I, II and III of Tabled re¬ 
spectively. Solid (dotted) arrows denote positions of xo (a:i) 
in Table [D Results of systems with different sizes Ly are al¬ 
most degenerate. 



yo/t 

Oijt 

Xq j Lx 

xxjLx 

Ax 

CASE I 

0.8 

0.0015 

0.458 

0.0833 

21 sites 

CASE II 

0.5 

0.0015 

0.375 

0.0833 

10 sites 

CASE III 

1.25 

0.005 

0.333 

0.0833 

4 sites 


TABLE L Three parameter setups we chose, GASEs I, II, and 
III. The last three columns xq/L^, x\ j and Ai denote rep¬ 
resentative positions in the incommensurate filling (IC) and 
the Mott regions, and the width of the IG region, respectively 
(see the text and Fig. EJ. 


Boundary state in T ^ 0 limit.— In order to un¬ 
derstand the boundary nature of the Bose Mott state, 
we will study particle densities, superfluidity, and cor¬ 
relation functions of the model m by QMC computa¬ 
tions. In Fig. [21 we first show the local density profile 
n(x) = {nx,y) and local helicity modulus Y{x), which 
is proportional to the local superfluid density |^, 
changing fj,Q and the curvature a at very low tempera¬ 
tures T {ks is set to be unity). We here show QMC 
results of three parameter settings {y,o,a) in Table [I] as 
representatives. In all CASEs I, H and HI, an incom¬ 
mensurate (IC) filling region (e.g., 0.21 < xjLx < 0.6 
in CASE I) appears between the filling-one Mott and 
the vacuum (empty) states. In the IC region, the lo¬ 
cal superfluid density takes a finite value. Size of the 
IC region and the superfluid density profile are almost 
irrespectively of the length Ly when Ly is sufficiently 
large. We note that finite-size effect along x direction can 
be ignored in the present parameter settings. Survival 
of a superfluid density sharply contrasts with the case 
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of the purely ID Bose system because the latter super- 
fluid density is known to disappear in the thermodynamic 


limit 331-33. We confirmed that a further narrower IC 


region still survives when a more large value of a is ap¬ 
plied. Furthermore, we observe an IC region when we 
apply other types of the chemical potentials with non¬ 
harmonic curvatures [si [S^: V{x) = and 

= ^2 + 0.2 exp(—a;/.^2)- These results indicate that an IC 
superfluid region between Mott and vacuum areas gen¬ 
erally appears and a special potential V(x) with an ex¬ 
tremely large curvature is necessary to remove the IC 
region. 
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FIG. 3: (Color online) Off-diagonal and density correlation 
functions along the y direction. Pairs [(a),(b)], [(c),(d)] and 
[(e),(f)] are the results of CASE I, II, and III at the IC po¬ 
sition X = xo, respectively. Pairs [(g),(h)] and [(i),(j)] are 
respectively the results of spatially uniform ID and 2D Bose 
Hubbard models at an IC filling. Black lines are guides to 
eyes. 


In Fig. [31 we present the equal-time one-particle (off- 
diagonal) correlation function C°{y)={b\.^ yhxQfi) and 
density one C‘^{y)={nxa,ynxo,o) - {nxo,y)^ at the bound¬ 
ary position X = xq where Ps{xq) takes the largest val¬ 
ues. As a comparison, we also show the QMC results 
of a purely ID Bose Hubbard model [Fig. |31[g) and (h)] 
and a spatially uniform 2D Bose Hubbard one [Fig. [31[i) 
and (j)] in an IC-filling case. The ID and 2D Bose sys¬ 
tems of Fig. [3] belong to TLL and Kosterlitz-Thouless 
(KT) phases, respectively. At the position xq, power- 
law decays along the y direction are observed in both 
off-diagonal and density correlations with long distances. 
Their critical exponents are evaluated by assuming the 
form C^{y)=const. x y~'^’‘, and they are summarized in 
Table m The emergence of algebraic decay is indepen¬ 
dent of the width of IC region. This clearly indicates 
that at least one gapless edge mode around the bulk 
Mott-insulating region always appears at very low tem¬ 
peratures. The algebraic decay of the density correlation 
is quite different from that of the KT phase in uniform 


CASE I 

CASE II CASE III 

ID 

Po 0.069(2) 

0.12(4) 

0.25(1) 

0.443(8) 

Pd 2.9(1) 

3.1(2) 

2.1(1) 

2.23(6) 


TABLE II: Estimated decay exponents for off-diagonal and 
density correlations, assuming C^{y) ~ const, x y~'^‘, where 
z = o (d) stands for the off-diagonal (density) correlation. 


2D Bose systems. In fact, Fig.[31[j) shows that the density 
correlation decays exponentially in the 2D case. In addi¬ 
tion, it is known 23 that two critical exponents satisfy 
Poyd = 1 in the purely ID TLL (see Table El, while the 
relation is clearly broken in the present boundary gap¬ 
less mode. These results conclude that the boundary IC 
region possesses intermediate properties between ID and 
2D Bose systems. 

Boundary state in finite temperatures.— Next, we dis¬ 
cuss the temperature dependence of the boundary IC 
states. In the present model as well as real experimen¬ 
tal systems, effects of finite size and spacial inhomogene¬ 
ity may spoil true phase-transition phenomena in the 
thermodynamic limit, but their residual things might 
still survive. Figure [3] shows the temperature depen¬ 
dence of local helicity modulus K(xo) at a typical IC 
position X = Xo in CASE 1. The figure shows that the 
Ly dependence of K(xo) becomes negligibly small below 
Tjt ^ 0.075, where off-diagonal and density correlation 
functions decay algebraically. This small Ly dependence 
implies the existence of a gapless KT-like phase around 
the bulk Mott state. 

To quantitatively determine a KT-transition-like tem¬ 
perature in our finite inhomogeneous system, we simply 
apply the standard finite-size analysis for the KT transi¬ 
tion in spatially uniform 2D systems, which is expected to 
be reliable if the width of the IC region is sufficiently large 
as in CASE 1. In the uniform system, Y approaches to 
2 A:bTxt/ 7I' at the KT transition temperature T = Tkt- 
For each finite system with size L, the KT transition 



FIG. 4: (color online) Local helicity modulus at a IC position 
X = xo in CASE 1. 

temperature T*{L) is known to be 

T-(L)=r,Hi^oo)(l + 5i^), (2) 

where C is a fitting parameter. Since we have the data 
of Y{Ly) for different sizes Ly in Fig. [31 the tempera- 
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FIG. 5: (a) Size dependence of characteristic temperature 
T*{Ly). (b) Temperature dependence of the exponent r^o- 


ture T*{Ly) can be estimated as the cross point between 
numerically determined Y{Ly) in Fig. 2] and the linear 
line Y{T) = 2kBT/Tr. In Fig. [H^a), we plot the eval¬ 
uated T*{Ly). In its inset, we determine the value of 
T*{Ly —oo) by combining T*{Ly) and the scaling re¬ 
lation ([2]) (T*{Ly —>• oo) corresponds to the KT transi¬ 
tion temperatnre Tkt in the case of uniform systems). 
The characteristic temperature T*(oo) is determined as 
T*{oo)/t ~ 0.11(1). 

As an alternating way to determine T*(oo), we can uti¬ 
lize correlation functions in our inhomogeneous system. 
It is well known that the critical exponent rjo of the off- 
diagonal correlator increases from zero and it becomes a 
quarter at the KT transition with the growth of tempera¬ 
ture. Let us simply apply this property to fix the KT-like 
temperature in our IC region. The inset of Fig. [Hb) is 
the T dependence of the exponent r]o. We see that rjo 
indeed crosses a quarter around T = T*{oo) which was 
estimated above from Y{Ly). 

We stress that the temperature T*(oo) is much lower 
than the true KT transition temperature Tkt in the uni¬ 
form 2D system. From the QMC calculation, we obtained 
TKx/t ^ 0.92(3) for U/t = 2Q and ii/t = —0.74, where 
the averaged particle number per site is almost same as 
that at the position x = xq in CASE I. This must be be¬ 
cause the development of off-diagonal correlation along 
the X direction is suppressed owing to the existence of 
Mott-insulating and vacuum regions. When the width of 
the IC region is small as in CASE III, it is hard to quan¬ 
titatively determine T* . However, even in snch a case, 
a KT-like power law in correlations appears at very low 
temperatures (see Eig. El). 

Structure factors.— Erom all the discussions above, we 
see that at least a gapless IC state always appears around 
the Bose Mott state if temperature is low enough. Fi¬ 
nally we discuss a experimental way to detect the gap¬ 
less edge mode. In cold-atomic systems, the momentum 
distribution of correlation functions [s^ can be observed 
in principle. For example, time-of-flight (TOF) method 
and light-scattering spectroscopy have been applied to 


observe them. In Fig. [6l we show the momentum-gy 
distribution of S°{x,qy) = C'°(a:, y)e“*^'^'' for 

the IC region at x = xq and for the bulk Mott region 
at a; = xi. Here, C°{x,y) = {hX^yb^fi). In the realistic 
experimental setup, the number of sites along the y-axis 
is less than ^100 sites and then we set Ly = 64 in all 
the panels of Fig. |6l As temperature decreases, the mo- 
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FIG. 6: (color online) (a)-(c) T dependence of structure fac¬ 
tors S°{x,qy) at Ly=64:. Gircles (triangles) are the results for 
the IC (Mott) region at x = xo (x = xi). (d) T dependence 
of the structure factor of the off-diagonal correlator S°{q) in 
a TLL state of ID Bose Hubbard model. 


mentum distribution at zero wave number = 0 of the 
IC region well develops, while that of the Mott region 
is suppressed for any wave number g^. The g^ = 0 peak 
reflects the development of superfluidity in the IC region. 
In Fig. E] (d), as a comparison to the IC region, we show 
the momentum distribution of a finite-size ID Bose Hub¬ 
bard model under an uniform chemical potential with 
the almost same filling as the IC boundary region. In 
the ID case, we also observe a py = 0 peak structure 
which is the contribution of TLL. From Fig. EJ we find 
that momentum distributions in both Mott and IC re¬ 
gions exhibit the similar T dependence to the finite-size 
ID system with the same filling. This is an another evi¬ 
dence for the existence of a gapless edge mode in the IC 
region and it also indicates the difficulty of distinguishing 
the IC gapless state and the ID TLL. 

Summary and discussions.— In conclusion, we have 
studied the edge state surrounding the 2D Bose Mott- 
insulating phase. From the QMC method, we have found 
that in the IC edge region, both off-diagonal and den¬ 
sity correlators show an algebraic decay and a snperflnid 
density appears below a characteristic temperature irre¬ 
spective of the width of the IC region. This ’’universal” 
gapless edge mode can be detected e.g., by observing a 
Qy = 0 peak of S°{x,qy). 

Our result naturally indicates that a similar gapless 
edge mode generally emerges in any kinds of 2D cold- 
atomic Bose and Fermi insulating states. Therefore, if 
we consider a topological insulating state in 2D cold- 
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atom systems, we can expect that there are at least two 
gapless modes of the edge state in the IC region and 
an intrinsic topologically-protected edge state as shown 
in Fig. B There might be a relevant coupling between 
these two edge states. Thus, when we discuss a way to 
detect topological edge mode in 2D cold-atom systems, 
we should generally consider effects of a non-topological 
(but universal) edge mode in the IC region around the 
bulk insulating area. 

(a) (b) 



FIG. 7: (color online) Spatial strnctures of a trapped Bose 
Mott state (a) and a trapped topological insulating state (b) 
in 2D cold atoms. 
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